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Without our ability to model and manipulate the band structure of semiconducting materials, the 
modern digital computer would be impractically large, hot, and expensive. In 8.06, we studied the 
effect of spatially periodic potentials on the spectrum of a charged particle in one dimension. We 
would like to understand how to extend these methods to model actual crystalline materials. Along 
the way, we will explore the construction of periodic potentials in three dimensions, and we use this 
framework to relate the single-particle Hamiltonian to the potential contribution from each atom. We 
then construct a crude model system analogous to the semiconductor silicon, and demonstrate the 
appearance of level splitting and band gaps as the strength of the potential is varied, in accordance 
with our intuition from the one-dimensional case. We discuss refinements of the model to include 
many-particle effects, and finally we show how a careful choice of the potential function leads to 
good agreement with the correct band diagram for silicon. 



I. INTRODUCTION 

Many world-changing results from solid state physics, like 
the development of the silicon transistor, take advantage 
of the band structure of crystalline materials to give en- 
gineers detailed control over material behavior. Because 
crystals are ultimately many-particle systems with many 
degrees of freedom, the models that are used to describe 
these materials make a wide range of approximations 
in order to obtain mathematically and computationally 
tractable results about band structure. The most detailed 
models can predict not only the number and spacing of 
bands and the density of available states in the material, 
quantities useful for understanding the behavior of elec- 
trons near the Fermi surface at nonzero temperature, but 
also the effects of mechanical strain, defects, and impuri- 
ties on the crystal's electronic properties. Semiconductor 
engineers fit these models to empirical data to obtain very 
nuanced control over the behavior of devices. 

In 8.06, we explored the behavior of single electrons 
in one-dimensional periodic potentials. This toy model 
allowed us to invoke several of the most important tech- 
niques for learning about band structure, including the 
nearly-free electron model, which explained how tiny in- 
teractions between otherwise free valence electrons can 
open up band gaps, and the tight-binding model, which 
showed that even closely-held core electrons undergo level 
splitting to form narrow bands. 

In this work, we will show how these tcchniqiics gen- 
eralize from one dimension to three dimensions, demon- 
strating the appearance of level splitting and band gaps 
as a periodic potential is gradually "switched on". We 
will also show how an approximation to the correct band 
structure of the semiconductor silicon can be obtained 
within this model by adjusting the potential function. 

In order to apply the nearly-free electron model, we will 
need a mathematical framework for the three-dimensional 
lattice. We can then take a periodic potential on this lat- 
tice over to the Fourier domain, recovering the reciprocal 
lattice. This will tell us how the Brillouin zone generalizes 
to three dimensions. After reviewing Bloch's theorem, we 



will put all of these pieces together to find and solve the 
Hamiltonian for a particle moving in a three-dimensional 
periodic potential, and consider briefly how this approach 
can be extended to the many-particle case. 

II. LATTICE FORMALISM 

A. Bravais Lattice 

Consider the potential energy of an electron located at a 
point X within an infinite grid of positively-charged ions, 
all carrying the same charge. We can write this as 

Vix) = J2U{x-Xi) (1) 

i 

where Xi is the position of the i^^ atom, and U{x) is 
the potential due to a single ion at the origin. In one 
dimension, we write the positions of the atoms as 

Xn = nax = na\ (2) 

for integer n, where o is the lattice spacing, and we have 
defined di = ax. The region of 1-D space bounded by 
(0, a) is a unit cell. 

Notice that there is no particular reason to require that 
U{x) be the potential due to a single ion. One could 
construct a 1-D crystal where small, identical clusters of 
atoms form a repeating pattern with period a, with each 
cluster containing several atoms spaced by less than a. 

In three dimensions, the situation is more interesting, 
since we have some freedom to define the directions along 
which the crystal is periodic (they don't necessarily have 
to be perpendicular). Even for oblique directions of peri- 
odicity, we can write 

Xn,m,e = nai + ma2 + ias (3) 

to give the positions of a regular grid of points. Families 
of points describable in these terms are called Bravais 
lattices. Common Bravais lattices include the cubic lattice 
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system (simple, body-centered, and face-centered) and 
the hexagonal lattice. The cubic lattice system is shown 
in Figure [T]^a) , (b), and (c), with circles indicating the 
locations of atoms in the simple, body-, and face-centered 
cubic configurations, respectively. You should take a 
moment to convince yourself that shifting the origin by 
one step along any of the indicated directions of periodicity 
«i,2,3 leaves the lattice unchanged. 

As in the 1-D case, it should be stressed that the poten- 
tial U {x) can describe a cluster containing more than one 
ion. For instance, diamond, silicon, germanium, and the 
allotrope a-tin adopt a crystal structure that positions 
two atoms at each point of a face-centered cubic (FCC) 
lattice. This structure, known as "diamond cubic", is 
shown in Figure [ij^d) . 



B. Bloch's Theorem 

If a Hamiltonian commutes with a set of unitary, mutually 
commuting translation operators Tj, then the eigenstates 
of H can be chosen to be (simultaneously) eigenstates 
of all the T^. This follows because hermitian and uni- 
tary operators are normal operators, and a set of normal 
operators commutes if and only if it is simultaneously diag- 
onalizable (an extention of our 8.05 result for commuting 
hermitian operators). Now any two translation operators 
of the form = exp(— • a/h) commute, since the indi- 
vidual components of p commute. Furthermore, since the 
translation operators are unitary, their eigenvalues have 
unit magnitude. 

Suppose that H can be written as the sum of a kinetic 
term /2m and a potential term which satisfies V{x) — 
V{x + a) for some set of a. Then [H, Tg] = 0, and we can 
find an eigenbasis for Ji which satisfies 

TM^e^'^m (4) 

for a in our set. 

We know that composing two translation operators Ts 
and Tg gives To maintain consistency with Q, 6s 

must be linear in a. We will take 

9s = -K-a (5) 

This gives 

(xifkiV') = e-'^-'im 

iPix-d) = e-"^-^ij{x) 

So we see that the eigenstates of H can be written in the 
spatial basis as 

^{x) = e^^-^ti(x) (6) 
where u{x) = u{x -f- a) for all a in our set. 



C. Reciprocal Lattice 

When we apply Bloch's theorem to the potential ([T]), we 
find that the wavefunction can be written in terms of 
a function u{x) which has the same periodicity as the 
lattice: 

u{x) — u{x + noil + Tiia2 + ^03) V n, m, ^ G Z (7) 

It is productive to ask whether the Fourier transform of 
such a function (that is, its representation in the momen- 
tum basis) has any helpful properties. Starting with 

u{x) = u{x + a) 

for some a, and writing each side in terms of the inverse 
transform of u{k), we find 

Collecting terms and taking the inverse transform gives 

0^(1- e'^-'^)u{k) 

from which it follows that 

Vfc, either a ■ k — 2ttN for some iV e Z, or u{k) — 0. 

This result shows that if a function is periodic in a, then 
its transform must be zero everywhere except where a ■ k 
is a multiple of 27r. This can be summarized by saying 
that periodic functions have discrete transforms. 

In three dimensions, we have three linearly independent 
directions of periodicity, and hence three discreteness 
constraints on the k where u{k) is allowed to be nonzero. 
Writing a - k as aJk, we can combine the three constraints 
into one matrix expression: 
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If we define 51,2.3 to be 27r times the three columns of 
the inverse matrix, then we have an expression for a new 
Bravais lattice: 

k = Ngi+Mg2 + Lg3 (8) 

Each k vector allowed by this lattice corresponds to a 
complex exponential e^^'^ in real space which is periodic 
in ai_2,3. This structure is called the reciprocal (or dual) 
lattice, and the space it occupies has many names: it is 
alternatively called reciprocal space, dual space, fc-space, 
transform space, or (loosely) momentum space. From 
now on, we will refer to the original Bravais lattice as the 
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FIG. 1: Cubic Bravais lattice system: (a) simple cubic, (b) body-centered cubic, (c) face-centered cubic. The vectors ai,2,3 
represent directions of periodicity. The unit cells contain 1, 2, and 4 points, respectively, (d) Diamond cubic variation, with 
additional atoms highlighted. The lines illustrate the tetrahedral symmetry about each atom. Atoms in (d) connected by a 
thick line are associated with the same lattice point. 



real space lattice. Curiously, if the real space lattice is 
face-centered cubic, the reciprocal lattice turns out to be 
body-centered cubic, and vice versa. 

The preceding argument established that the Fourier 
transform of a periodic function is nonzero only on the 
reciprocal lattice. If we write the potential ([T]) as a con- 
volution product, 

V{x) = U (x) * S{x — ndi — ma2 — ia^) 



we can take the Fourier transform of the two factors 
separately to obtain 



V{k) = U{k) J2 



-ik-(nai-\-ma2-\-ia3) 



When we carry out the summation, we get a new sum 
over the reciprocal lattice, which we can again write in 
matrix form: 



N.M.L 



(27r)3 6{k ■ ai-2TrN)S{k ■ a2-2iTM)5{k ■ a3~2nL) 
= (27r)3 ^ 5 
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Finally, we can show by u-substitution that this is just 
{2nr 



det \a1a2a3 



Y,S{k-G) 



where is over the reciprocal lattice. So 

v{k) u{k) f '^(^ - ^) 

det 1010203! t;^ 



[2^ 



J2UGS{k-G) 



(9) 



G 



where we have defined Ug = U{G) and D, = det |aia2a3| 



space. But we just observed that every e^^'^ with G on the 
reciprocal lattice is periodic in real space. Thus if we write 
ip'{x) — e'(''~'^)'^(e''^'^u(a;)) — e^'^' u' {x) , we obtain an 
equally valid expression for the exact same wavefunction. 
In this sense, k is not really a free parameter of the system, 
but can be restricted without loss of generality to lie in a 
small region of the reciprocal space. Conventionally, we 
choose this region to contain all vectors k which are closer 
to G = than to any other reciprocal lattice point. If we 
start with a k closer to some other reciprocal lattice point 
G ^ 0, we can subtract G from it to get an equivalent 
point closest to 0. The region of space closer to than 
to any other reciprocal lattice point is called the first 
Brillouin zone. In one dimension it was given by the 
interval (— -, -), since the boundaries between regions 
in one dimension are simply points. In three dimensions, 
these boundaries are planes, and the first Brillouin zone 
has a rather beautiful shape formed by the intersections 
of boundary planes. The first Brillouin zone for an FCC 
real space lattice is shown in Figure [2j This polyhedron 
bounds the region of reciprocal space where we will focus 
our attention. 

The job of a band structure diagram is to summarize 
all the interesting things going on in fc-space. Since it 
is difficult to capture all the features of a volume on a 
fiat piece of paper, these diagrams typically take a little 
tour around the first Brillouin zone, following a (some- 
what) standardized trajectory. You will notice that a 
number of particularly symmetric points in the figure 
have been named with capital Roman and Greek letters. 
These letters will appear along the horizontal axis of band 
structure diagrams, and it is understood that motion 
along the horizontal axis corresponds to progress along a 
trajectory between these points of symmetry in fc-space. 
For instance, the segment of the graph labeled T-X cor- 
responds to the straight line in fc-space from /c = to 
k — {2iT/a)x. 



D. First Brillioun Zone 



III. NEARLY-FREE ELECTRON MODEL 



Our study of Bloch's theorem led us to wavefunctions 
of the form ^{x) = e*'''^M(x), with u{x) periodic in real 



We would like to write down and solve the Hamiltonian 
for a single particle moving in a periodic potential. While 
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can write the eigenstates \'ip) as 



|Wk) 



(13) 



where Uft(x) has the same periodicity as the real space 
lattice, and the Schrodinger equation is written as 



1 



n,s^ —{p + MY + V{x) 
2m 



(14) 
(15) 



Using the plane wave basis for we find matrix ele- 
ments 



FIG. 2: First Brillouin zone for a face-centered cubic real space 
lattice. Represents the volume in fc-space closer to than to 
any other reciprocal lattice vector, and includes all the Bloch 
wave vectors k necessary to describe the band structure of the 
crystal. 



it is not true in general that the dynamics of a many- 
electron system can be described in terms of one set of 
single-particle states, single-particle wavefunctions are 
nevertheless a useful tool for approximating the correct 
many-particle system. 

Working in the plane wave basis, we will find an expres- 
sion for the matrix elements of the Bloch Hamiltonian in 
terms of the Fourier transform of the potential. We ne- 
glect electron-electron interactions. Later on, we will show 
how to interpret this approximation as an application of 
the variational principle. 



A. Basis 

Our orthogonal basis and normalization convention are 
{x\x'} = 6{x-x') 

d^xe'^-^\x) (10) 

d^x J d3f'(f|e-''*^-V'='-^'|x') 
{2TrfS{k-k') (11) 



\k) 
{k\k') 



B. Hamiltonian 

We begin with the Hamiltonian for a particle in an arbi- 
trary potential. 



(12) 



Using the fact that 



for any translation a in 



(fclHeff |fc') = —ik + ^f{k\k') + {k\ViS)\k') (16) 

Zm 

We can write V{x) in terms of V{x) with the help of (|9|: 
1 



^ (2^ y '^^^^^^^ 

= ^^Ug I d^kS{k ~ G)e'^-'^ 



Using the sifting property of the delta function. 



1 



iG-x 



We can now find matrix elements for our potential using 
the definition of |fc): 



{k\vCx)\k') = -J2UGl d^x I dV 



G 



(x\e-'^-'^e'^-^e'^'-^'\x') 



G 



J2Ug d'^xe'" 



x-(G~k+k') 



Taking the integral, we have 

{k\Vmk') = h2^fY.UG5{G-k + k') 



(17) 



Computations in a continuous basis like \k) are problem- 
atic because they are difficult to approximate numerically. 
We would much rather work in a vector space with a dis- 
crete basis, so that our linear operators become ordinary 
matrices. Fortunately, the structure of pT] ) permits us to 
divide the continuous basis into a set of discrete subspaces, 
with the assurance that there will be no coupling between 
subspaces. In matrix terminology, our Hamiltonian is 
block diagonal. Using this fact, we can easily restrict the 
action of the Hamiltonian to a given subspace and solve 
numerically. 



the real space lattice and applying Bloch's theorem, we 



From (17), we can see that the Hamiltonian (15) cou- 
ples states and \k') if they differ by a reciprocal lattice 
vector. Otherwise, the off-diagonal matrix elements van- 
ish. Thus, to every k in the first Brillouin zone there 
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corresponds a family of states \k + G) which are (po- 
tentially) coupled by the potential. Without losing any 
physics, we can restrict the Hamiltonian to act in one 
of these discrete orthonormal subspaces by constructing 
(and normalizing) matrix elements 



{k + G\Vix)\k + G') 
{k + G'\k + G'){k + G\k + G) 



V(27r)35(0)(27r)3,5(0) 
1 

n 



Ug-g' 



(18) 



We can build one of these matrices for each k we are 
interested in, add h?{k + G)^/2m to the diagonal, and 
we will have the full Hamiltonian for the corresponding 
subspace. 

Now that we have matrix elements, why don't we apply 
perturbation theory? If the potential is small compared 
to h?(k + k)^/2to, we can view it as a perturbation on 
the Hamiltonian for a free particle. In order to obtain 
corrections to the free particle energies and wavefunctions, 
we must first look for the "good" states, that is, the basis 
in which the perturbation is diagonal within every de- 
generate subspace of the unperturbed Hamiltonian. We 
know that the degenerate subspaces of the unperturbed 
(free) Hamiltonian correspond to spheres of constant |fc| 
in momentum space. The question is whether the matrix 
element of V{x) between two states on the same sphere 
is nonzero. Since a sphere has chords of every length 
in (0, 2r], we can see by (181 that any vector G in the 



reciprocal lattice which is short enough to fit inside the 
sphere (that is, with < |G| < 2\k\) couples together 
degenerate states. Certain special potentials, for instance 
those which are periodic in x but constant in y and z, per- 
mit us to choose a sine and cosine basis that diagonalizes 
V{x). For more general potentials, finding a good basis 
is bothersome. At this point, the easiest thing to do is to 
throw the problem at a computer. We have to make a few 
more decisions - most importantly, what our potential 
function U{x) is, and how many states \k + G) to include 
in the truncated (finite) matrix form of JieS ~ but given 
these, finding the single-particle energies and wavefunc- 
tions is a simple matter of computing each matrix element 
numerically using ( |18[ ), and then invoking a subroutine 
to find the eigenvalues and eigenvectors of this matrix. 
Fast and accurate subroutines to perform this function 
are present in every standard linear algebra system. In 
the popular LAPACK library (Linear Algebra PACKage), 
for instance, there is a subroutine called zheevd that does 
what we want. To generate a band structure diagram, we 
will need to compute and plot these eigenvalues for each 
k we encounter on our tour of the first Brillouin zone. 



IV. MODEL SYSTEM 

Using the framework we have developed, we can now 
choose a particular system and plot some band diagrams 



to verify our work and improve our intuition into the 
phenomenology of band structure. 



A. Potential 

As it turns out, a great deal of research [1-5 has gone 
into the choice of the potential function, since much of 
the actual many-particle physics can be incorporated by a 
suitable choice of U{x). For instance, once we obtain wave- 
functions for a crystal filled with non-interacting electrons, 
we can calculate the effect of e~-e~ repulsion (classical) 
and correlation (the quantum "exchange force"). Treating 
these as an adjustment to the potential, the system can 
re-solved to obtain a set of wavefunctions that are sensi- 
tive to many-particle effects. These can in turn be used to 
make a better estimate of repulsion and correlation effects. 
The iterations stop when the wavefunction-dependent po- 
tential and the potential-dependent wavefunctions are 
mutually consistent. This is the Hartree-Fock method. 
The approximation here lies in the assumption that the 
many-particle ground state can be written as an anti- 
symmetrized product of single-particle wavefunctions. In 
general, it is a linear combination of many such prod- 
ucts. This makes Hartree-Fock a variational technique, 
since we are effectively minimizing the energy of a trial 
multi-particle wavefunction. The minimization is over all 
possible sets of single-particle states. Since the expecta- 
tion value of the Hamiltonian on any state is an upper 
bound on the ground state, the Hartree-Fock method 
gives an upper bound on the energy of the crystal system. 

For simplicity, we will use a Coulomb potential of the 
form 



U{x) = 



\r\ 



(19) 



Since our goal is to model the diamond-cubic silicon crys- 
tal, we will have to add another Coulombic potential at 
an offset to represent the two-ions-per-lattice-site struc- 
ture. Even this gross simplification (we haven't consid- 
ered e~-e~ interactions, spin, etc.) is sufficient to capture 
some of the flavor of the correct band diagram. In any 
event, we'll need the Fourier transform of this potential, 
so we'd better go ahead and compute it. 

As it turns out, the Fourier integral for the function 
— Zasq^/r doesn't converge. The transform of the related 
"Yukawa" potential —Zcsq^e~^^^/r is 



u{k) = -z,m' 



fc2 



(20) 



In the limit as to —j' 0, we recover the transform of the 
Coulomb potential everywhere except at fc = 0. 

,47r 



Uik) = ~Z,m' 



(21) 



Since Uq only appears on the diagonal of the Hamiltonian, 
ignoring the divergence at /c = corresponds to neglecting 
a constant energy shift. 
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B. Results 

Computed band structure diagrams for the diamond cubic 
crystal (Coulombic potential) are presented in Figure [s] 
for four different values of Zgg. The labels along the 
horizontal axis correspond to fc-space points in the first 
Brillouin zone (Figure [2]) . Note especially that 

• For a free electron, we expect E = /2m cx fc^. In 
fact, when Z^g = 0, we find parabolic dispersion 
around the point fc = F = 0. 

• Along the path L-F, we recover much the same 
"folded parabola" structure we saw in one dimension, 
with branches of the band diagram reflecting off the 
wall of the Brillouin zone (in this case, at L). 

• For small Z^g, the effects of the perturbation 
should be most visible near crossings and degen- 
erate branches of the band diagram, where level 
repulsion is strongest. This effect is observed as 
expected. Note that even though the Coulomb po- 
tential gives nonzero Fourier components for every 
nonzero k, not all degeneracies are broken. Part of 
the reason is due to zeros imposed by the two-atom- 
per-FCC-site diamond cubic structure. 

• For a tightly bound electron (high Z^s), we expect 
to find more or less isolated bands corresponding to 
the orbitals of isolated atoms. Since band spreading 
corresponds to tunneling in the tight binding picture, 
these bands describe localized states. 

In Figure |4j we show how this simple model can incor- 
porate the more complicated situation of many-particle 
interactions in an ad-hoc way. By selecting a few matrix 



elements essentially by hand, a very close correspondence 
can be obtained with results from more advanced litera- 
ture 1; . In the literature, these matrix elements would be 
calculated, if not from first principles, then from fitting 
with experimental data such as x-ray reflectivity curves. 



V. DISCUSSION 

As we have seen, it is possible to get a lot of mileage 
out of the single-particle system. By adopting more so- 
phisticated approaches to the potential due to e~-e~ 
interactions, [T] and many others have obtained models 
for band structure in real semiconductors with good accu- 
racy and substantial predictive power. The single-particle 
framework we have developed is simplistic insofar as it 
ignores the effects of spin-orbit coupling, electron corre- 
lation, strains, electric and magnetic fields, impurities, 
and finite-temperature effects. It is nonetheless powerful 
enough to make contact with the more advanced litera- 
ture when supplied with an ad-hoc choice of potential, 
and it remains beautiful in its simplicity. 
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FIG. 3: Level splitting and band gap formation in a diamond cubic crystal as the strength of the Coulombic potential is varied 
upwards from zero. From top left to bottom right, Z^b = 0,0.25,0.5,2.0. Energy gaps are shown in gray. Letters along the 
horizontal axis represent points in the first Brillouin zone (Figure [2|. 
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(a) Selected matrix elements 
Vfc k+G foi' 3- non-Coulombic 
potential, where \G\ = nvr/a. 

Units are electron volts. 
Remaining matrix elements 
are unchanged. 
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(b) Computed diamond cubic crystal band 
structure using the non-Coulombic potential 
described in|(a)| 
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(c) Band structure for Si. Solid and dashed 
lines correspond to different techniques for 
incorporating e~-e~ interations. Reproduced 
with permission from 1 . Copyright 1976 by 
The American Physical Society. 



FIG. 4: By modifying the potential, good agreement can be obtained with the band structure for Si shown in Figure 1 of ^Ij. 



